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Worm Algorithm for Continuous-Space Path Integral Monte Carlo Simulations 
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We present a new approach to path integral Monte Carlo (PIMC) simulations based on the worm 
algorithm, originally developed for lattice models and extended here to continuous-space many- 
body systems. The scheme allows for efficient computation of thermodynamic properties, including 
winding numbers and off-diagonal correlations, for systems of much greater size than that accessible 
to conventional PIMC. As an illustrative application of the method, we simulate the superfluid 
transition of ''He in two dimensions. 
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Over the past two decades, PIMC simulations have 
played a major role in the theoretical investigation of 
quantum many-body systems, not only by providing re- 
liable quantitative results, but also by shaping our cur- 
rent conceptual understanding, e.g., of the relationship 
between superfluidity and Bose condensation. At least 
for Bose systems, PIMC is the only presently known 
method capable of furnishing in principle exact numerical 
estimates of physical quantities, including the superfluid 
density, and the condensate fractional . 

As PIMC is currently the most realistic option to in- 
vestigate ever more complex quantum many-body sys- 
tems, the issue arises of overcoming its present limita- 
tions. Aside from the notorious sign problem, the main 
bottleneck of the current PIMC technology is inarguably 
the maximum system size (i.e., number N of particles) for 
which accurate estimates can be obtained, in a reasonable 
amount of computer time. Almost two decades since the 
pioneering work of Pollock and CeperleyQ, who simu- 
lated the superfluid transition in bulk liquid "^He on a sys- 
tem of iV=64 atoms, no further advance has been made, 
despite a hundredfold increase in computer speedQ- 

For such quantities as energy and diagonal correla- 
tions, one can often approach the thermodynamic limit 
(iV —^ oo), by studying systems comprising as few as ^ 
30 particles. However, accurate predictions of superfluid 
properties of liquids (and solids [sjl) require that the 
superfluid and condensate fractions, ps and no, be com- 
puted for large systems of significantly different sizes. In 
conventional PIMC, ps is obtained by means of the so- 
called winding number estimator)^, which can only take 
on a nonzero value if long permutation cycles of iden- 
tical particles occur in the system. Because the sam- 
pling frequency for such cycles decreases exponentially 
with N, ensuring the ergodicity of the algorithm becomes 
problematic 'F]. 

This hurdle seems difficult to conquer within any 
scheme formulated in the canonical ensemble, in which 
the winding number is "topologically locked" in the 
N ^ oo limit 0]. On the other hand, the same hur- 



dle has been completely overcome in quantum Monte 
Carlo simulations of lattice models. A lattice Path Inte- 
gral scheme based on an alternative sampling approach, 
known as worm algorithm (WA) 6], allows for efficient 
calculations of winding numbers and one-particle Green 
function G, for systems of as many as ~ 10^ particles0- 
A fundamental aspect of the WA is that it operates in an 
extended configurational space, containing both closed 
world-line configurations (henceforth referred to as Z- 
or diagonal configurations), contributing to the partition 
function Z , as well as configurations containing one open 
line (worm). The latter configurations contribute to the 
one-particle Green function; below, they are referred to 
as G- (or, off-diagonal) configurations. All topologically 
non-trivial modifications of world lines occur in the off- 
diagonal configurational space, where there are no con- 
straints; when the sampling process generates a diagonal 
configuration, the number of particles and the winding 
number are updated. 

In this Letter, we describe the extension of the WA 
to the PIMC simulation of quantum many-body systems 
in continuous space. Our novel PIMC implementation, 
while based on the same theoretical underpinnings dif- 
fers fundamentally from the "canonical" one0, both in 
the configuration space structure, as well as in the sam- 
pling method. Since the number of continuous configura- 
tion variables is no longer conserved, the new scheme nec- 
essarily belongs to the generic domain of diagrammatic 
Monte Carlo methods As an illustrative application 
of this method, we simulate the superfluid transition in 
liquid ■*IIe in two dimensions (2D), for systems with up 
to N—2500 particles, i.e. two orders of magnitude larger 
than in the most recent PIMC studv[lo|. In particular, 
we observe a dramatic speed-up in convergence of ps and 
G. 

We begin by reviewing conventional PIMC. One ob- 
tains averages of physical quantities (at a temperature 
T) over a set of many-particle configurations {R}, statis- 
tically sampled from a probability density proportional 
to p{R,R,f3) = {R\e-l^"\R), where (3 = l/T (we set 
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fcB=l) and H is the system Hamiltonian. The goal 
is achieved by samphng discrete many-particle paths 
X = (Ri, R2, ■ ■ ■ , Rp), periodic in the imaginary time 
interval /3 — Pt, from the probability density 



p{X) = e~^wn'°°(^J'^J+i 

tN 
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where po{Rj , Rj+i,T) = Hi^i Po(ry , r^j+i, r) is a prod- 
uct of N free-particle propagators, whereas U incorpo- 
rates correlations, both in space and in imaginary time, 
arising from interactions among particles. U is chosen so 
that, in the r — > limit, the distribution of configura- 
tions R visited by paths X reproduce i?, (3). Several 
choices are possible for [/, but our algorithm does not 
depend on its particular form. 

Eq. J^l implies the following configuration space struc- 
ture: N single-particle paths (world lines), labeled i = 
1,2, . . . , N , propagating in the discretizcd imaginary time 
t from ii = to ip+i = p. Each world line is formed by 
P successively linked "beads" labeled by the number of 
the corresponding time slices j = 1,2, . . . , P. The j-th 
bead of the i-th world line is positioned at r,y. The [3- 
periodicity implies that the (P+ l)-th bead of each world 
line coincides with the first bead of either the same, or 
another world line. 

The set of paths {Xi} is sampled by a Metropolis ran- 
dom walk through configuration space. In order to gen- 
erate from the current path Xi, a local space-time 
modification of Xi is proposed. The new path X* is then 
either accepted, Xi+i = X*, or rejected, Xi+i = Xi, 
based on according to the standard procedure 

The basic update Xi — > X* consists of deforming one 
or more randomly selected world lines, over a number 
1< m < P of successive links. In order to incorporate 
effects of quantum statistics, it is also crucial to allow 
groups of 1 < n < A'^ world lines to exchange: A mod- 
ified portion of a world line in the group will connect, 
m links later, to a different world line, among the n se- 
lected. Note that the number of world lines along X is 
always equal to N, in this scheme. Updates with arbi- 
trary exchange cycles ensure ergodicity of the algorithm. 
Typically, however, the acceptance rate for permutations 
is frustratingly low, particularly in the presence of repul- 
sive inter-particle potentials (e.g., in condensed helium), 
rendering the calculation inefficient, and impractical for 
large N. 

The WA described here has the same starting point, 
namely Eq. JQ), but with a crucial generalization of 
the configuration space, which now includes both the 
above-mentioned diagonal and off-diagonal paths, the 
latter corresponding to the representation (analogous 
to Eq. QJ) of the one-particle Matsubara Green func- 
tion G{r,t). Each off-diagonal configuration contains a 
worm, that is, a world line (on a /3-cylinder) with two 
ends — the "head" and the "tail" — corresponding to the 



Green function annihilation and creation operators, re- 
spectively. The two special beads at the open world line 
ends are named (for historical reasons) Ira (J) and Masha 
[M). Configurations in which 2 and M are located in 
space-time at points {vx,tx) and (r^,i^) contribute to 
G{yx — J^M,ti ^ ^m) with the weight defined in accor- 
dance with the generalized Eq. 

The sampling of paths {Xi } is implemented in WA ex- 
clusively through a set of simple, local updates evolving 
X (or, Ai) in space-time. The particle number becomes 
configuration- and time-dependent (there is one less par- 
ticle between 2 and A^, than in the rest of the path). In 
other words, the WA opens up the possibility to work in 
the grand canonical ensemble, with the chemical poten- 
tial /i being an input Darameter|l2j. 

Next, we describe the set of ergodic local updates 
which sample our extended configuration space, switch- 
ing between the Z- and G-sectors. Updates which change 
the number of continuous variables in X, are arranged in 
complementary pairs, satisfying detailed balance. Gen- 
eral principles of balancing complimentary pairs can be 
found in Ref. We have three pairs: Open/Close, In- 
sert/Remove, and Advance/Recede. Only the Swap up- 
date in the list below does not fall in this category, be- 
cause it preserves the number of variables, i.e., it is self- 
complementary. Naturally, all known standard tricks can 
be used, in order to enhance performance. 

(la) Open. This update is only possible if the config- 
uration is diagonal. 13J A world line (say the i-th) and a 
bead (say the j-th) are selected at random. A random 
number (rn—1) of beads, namely j + l, j + 2, . . ., j+m—1 
are removed, so that a worm appears with I at (Yij,tj) 
and at (ri.j+„j, tj+„j). Hence, the difference between 
the proposed new path, X* , and the previous one, X , is 
that instead of the i-ih world line there are now two new 
ones: The z-th world line (we retain the same label) now 
ends at the j-th bead (T) and the iQ-th. world line (we 
introduce a new label) corresponds to the piece of the 
original z-th world line starting from the (j -|- m)-th bead 
{M.)- The acceptance probability for this update is 



Pop = min<^ 1, 
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AC/ — /j,mr 



(2) 



Po(r.y,rij+„,mT) 

where AU ^ U{X) - U{X*), Nx is the number of world 
lines (particles) in the diagonal configuration X , and an 
arbitrary constant C controls the relative statistics of Z- 
and G-sectors. The number fh < P defines the interval 
for m: m € [l,Tn\. In practice fh is adjusted to ensure 
the desired acceptance rate. Due to the /3-periodicity, 
without loss of generality we assume that whenever the 
situation j + m > P occurs, the enumeration is shifted 
in such a way that j + m < P, and no ambiguity oc- 
curs. These definitions are common to all other moves 
described below, in which m and rh enter. 

(lb) Close. This update is only possible if the config- 
uration is off-diagonal. 13j Let 2 be the j-th bead of the 
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i-th worldline and Ai be the {j + m)-th bead of the iQ-th 
worldhne. If to > to, the move is rejected. If m < to, one 
proposes to generate a piece of world line connecting X 
to Ai, thereby rendering the configuration diagonal. The 
corresponding spatial positions of new (to — 1) beads, 
r^j+i, . . . , rij4_m_i, are sampled from the product of 
TO free-particle propagators IIJLi Po(i"i j+^-i, rjj+y, r). 
The probability to accept the move is 

1^ CNx*^m J 

If the move is accepted, the label is removed. 

(2a) Insert. The other way to create an off-diagonal 
configuration from a diagonal one is to seed a new to- 
link long open world line in vacuum. The number of 
links TO < TO and the position of M in space-time are 
selected at random. The spatial positions of the other to 
beads are generated from the product of to free-particle 
propagators. The move is accepted with probability 
Pin = min{l, CVP rhe^'^+''"'''}, where V is the sys- 
tem volume. 

(2b) Remove. The removal of the worm, i.e., the world 
line connecting A4 to 2, is attempted if its length (in the 
/3-periodic sense) is to < to. (If m > ffi, the proposal is 
rejected Ts^.) The acceptance probability for the move 
is Fj-n, — min{l, e^^~''™'^/CyF m}. We are in position 
now to select C. A natural choice would be C = l/VPfh, 
so that the probability to Open a worm of zero length is 
N/V = G(0,-0). 

(3a) Advance. This move advances X a random num- 
ber TO of slices forward in time. It is similar to Insert 
update in implementation. The acceptance probability 



is -Pad = min{l, 



}. Note that it is possible for X 



to advance past A4. 

(3b) Recede. Now X moves backwards in time (in the (3- 
periodic sense) by erasing m consecutive links, the num- 
ber 1 < TO < TO is selected at random. The acceptance 
rate is Pro — min{l, e^^"'^™'^}. If to turns out to be 
equal or larger than the total number of links between X 
and M. along the world line connecting them, the update 
is rejected 13j. 

(4) Swap. Let X be positioned on the z-th world line at 
the j-th time slice. (See Fig. n) Consider all the world 
lines intersecting the (j -I- TO)-th (in the /3-periodic sense) 
time slice and select one of them (labeled below with k) 
with the probability Tk = po(jCij ,'rk.j+m,iT^T)/'Si where 



(4) 



is the normalization factor (if the selected world line 
contains M at j'-th time slice, such that j' G + 
to], the move is rejected). A set of random positions 
rij+i, . . . , ri.j4_m_i is then generated as in the Close 
move, whereas the beads r^j+i, . . . , rfcj_|_.,„_i are all 
erased. X is shifted to r^j, while the world line i re- 
connects with the rest of the world line k, which implies 



re-labeling, as illustrated in Fig.^ The move is accepted 
with probability Pgw — min{l, e'^'^Si/Efc}. 

The Swap move generates all possible many-body per- 
mutations through a chain of local single-particle up- 
dates. Since no two particles need be brought within 
a distance of the order of the potential hard core, it en- 
joys a high acceptance rate, similar to that for the Ad- 
vance/Recede procedures. It must be emphasized that 
in our algorithm, unlike in conventional PIMC, arbitrary 
permutations of identical particles, as well as macroscopic 
exchange cycles appear automatically, if the physical con- 
ditions warrant them. This is because the statistics of the 
relative positions for the worm ends is given exactly by 
the Green function G{r,t). 
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FIG. 1; (Color online). Schematic illustration of Swap move 
described in the text, (a): before the move, (h): after the 
move. 

As an illustrative application of the new method, we 
present here simulation results of the superfluid transi- 
tion in 2D helium. Specifically, we have repeated the 
study first carried out in Ref. |l4| . using the same in- 
teratomic potential 15] and at the same 2D density 
p = 0.0432 A^^, but on systems with a number of parti- 
cles up to hundred times larger. We used an approxima- 
tion accurate up to ^^] for the high-temperature den- 
sity matrix [which determines the structure of the func- 
tion U(X)\, and extrapolated the results to the r ^ 
limit. For a given choice of r, the statistical error of ps 
is comparable to that of the kinetic energy. 

Fig. El shows our results for the superfluid fraction 
ps{T), obtained on systems comprising different num- 
bers N of atoms. Using the procedure illustrated in 
Ref. based on Kosterlitz-Thouless theory [l^, we 

have obtained numerical fits to our data, in the critical 
temperature range 0.65 K < T < 0.8 K. Our estimates 
for the values of the fitting parameter are d — 8.8 ± 0.5 
A for the vortex core diameter, and E—2.18 ± 0.04 K for 
the vortex energy, which lead to an estimate for the crit- 
ical temperature rc=0.653±0.010 K, significantly differ- 
ent from the previous result, 0.72±0.02 K, deduced from 
the iV = 25 dataQ. 

As mentioned above, our method gives easy access 
to the imaginary-time one-particle Green function, and 
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FIG. 2: Superfluid fraction ps{T) computed for 2D ''He on 
systems with different numbers TV of ''He atoms. The system 
density is p = 0.0432 A~'^. Dashed lines represent fits to 
the numerical data (in the critical region) obtained using the 
procedure illustrated in Ref . [l^ . The leftmost dashed line is 
the extrapolation to the infinite system. Open squares show 
results obtained in Ref. |3] for the same system, with N=2!3. 
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FIG. 3: (Color online) . One-particle density matrix computed 
for 2D *He at a density p = 0.0432 for a system of 

200 atoms, at r=0.675 K (upper curve) and r=I.O K (lower 
curve). Statistical errors on the curves are very small, and 
not shown for clarity. In the inset we present data (on a log- 
log scale) for the N=2500 system at r=0.675 K, with clear 
signatures of the Kosterlitz-Thouless behavior in the vicinity 
of the critical point. 



therefore to the one-body density matrix. For 2D hehum, 
this quantity is expected to decay to zero at all temper- 
atures, following a slow power-law behavior for T < Tc- 
Typical results obtained in this study, for systems with 



iV=200 and N = 2500 are shown in Fig. 01 

In conclusion, wc have implemented a novel procedure 
to perform large-scale PIMC simulations. Our scheme ex- 
tends to continuous space the worm algorithm previously 
developed for lattice systems, and affords efficient com- 
putations of thermodynamic properties, including the su- 
perfluid density and the single-particle Green function, 
for system of significantly larger size than accessible to 
the existing PIMC technology. 
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